#ifndef ARITHMETIC_H
#define ARITHMETIC_H

#include <boost/numeric/interval.hpp>
#include <fstream>
#include <vector>

using namespace boost::numeric;
typedef interval<double> interval_d;

class RandomInput
{
	int seed;
	static const int a = 1664525;
	static const int c = 1013904223;
public:
	RandomInput(int seed)
		:seed(seed)
	{ }

	// returns random integer in range 0 .. 2^31-1
	int get_int()
	{
		return (seed = seed * a + c) & 0x7fffffff;
	}

	double get_double()
	{
		const double lo = -1000;
		const double hi = 1000;
		const double scale = 1. / (1 << 30);
		return (lo + (hi - lo) * scale * (get_int() & 0x3fffffff));
	}
};


struct Point
{
	double x;
	double y;
	Point() { };
	Point (double x, double y): x(x), y(y) { }
	Point(RandomInput& rand):  x(rand.get_double()), y(rand.get_double()) { }
};

std::istream& operator>>(std::istream& in, Point& a)
{
	in >> a.x >> a.y;
	return in;
}

//adaptive arithmetic
std::pair<double, double> sum(double a, double b)
{
	double res = a + b, roundoff;
	double bv = res - a, av = res - bv;
	double br = b - bv, ar = a - av;
	roundoff = ar + br;

	return std::make_pair<double, double>(res, roundoff);
};

std::pair<double, double> split(double a, int s)
{
	double c = ((1LL << s) + 1LL) * a, ab = c - a;
	double hi = c - ab, lo = a - hi;

	return std::make_pair<double, double>(hi, lo);
};

std::pair<double, double> multiply(double a, double b)
{
	double res = a * b, 
		roundoff;

	int s = std::numeric_limits<double>:: digits / 2
		+ std::numeric_limits<double>:: digits % 2 ;

	double a_hi, a_lo, b_hi, b_lo;

	a_hi = split(a, s).first, a_lo = split(a, s).second;
	b_hi = split(b, s).first, b_lo = split(b, s).second;

	double e1 = res - a_hi * b_hi, 
		e2 = e1 - a_lo * b_hi, 
		e3 = e2 - a_hi * b_lo;

	roundoff = a_lo * b_lo - e3;
	return std::make_pair<double, double>(res, roundoff);
};

void grow_expansion(std::vector<double> &e, double &b)
{
	double q = b;
	for(int i = 0; i < e.size(); i++)
	{
		std::pair<double, double> s = sum(q, e[i]);
		e[i] = s.second;
		q = s.first;
	}
	e.push_back(q);
};

std::vector<double> expand_sum(const std::vector<double> &e, const std::vector<double> &f)
{
	std::vector<double> res = e;
	for(int i = 0; i < f.size(); i++)
	{
		double q = f[i];
		for(int j = i; j < res.size(); j++)
		{
			std::pair<double, double> s = sum(q, res[j]);
			res[j] = s.second, q = s.first;
		}
		res.push_back(q);
	}
	return res;
};

//left turn & bounding box

bool bounding_box(double a, double b, double c, double d)
{
	if(a > b)
		std::swap(a, b);

	if (c > d)
		std::swap(c, d);

	if(a > d || c > b)
		return false;
	return true;

};

int left_turn(Point a, Point b, Point c)
{
	double eps_m = std::numeric_limits<double>::epsilon();
	double eps = 4 * eps_m * std::abs((b.x - a.x) * (c.y - a.y)) + std::abs((b.y - a.y) * (c.x - a.x));
	double s = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);

	std::vector <std::pair<double, double> > mult;
	std::vector<double> e1, e2, e3;

	if(s > eps)
		return 1;
	if(s < -eps)
		return -1;
	
	interval_d inter_res = (interval_d(b.x) - interval_d(a.x)) * (interval_d(c.y) - interval_d(a.y)) -
		(interval_d(b.y) - interval_d(a.y)) * (interval_d(c.x) - interval_d(a.x));

	if(inter_res.lower() > 0)
		return 1;
	if(inter_res.upper() < 0)
		return -1;

	mult.push_back(multiply(b.x, c.y));
	mult.push_back(multiply(b.x, -a.y));
	mult.push_back(multiply(-a.x, c.y));
	mult.push_back(multiply(-b.y, c.x));
	mult.push_back(multiply(b.y, a.x));
	mult.push_back(multiply(a.y, c.x));

	std::vector<double> exp;
	for(int i = 0; i < 6; i++)
	{
		grow_expansion(exp, mult[i].first);
		grow_expansion(exp, mult[i].second);
	}

	for(int i = 11; i >= 0; i--)
	{
		if(exp[i] > 0)
			return 1;
		if(exp[i] < 0)
			return -1;
	}
	return 0;   

};

//intersection function
bool is_intersection(Point a, Point b, Point c, Point d)
{
  return bounding_box(a.x, b.x, c.x, d.x) && bounding_box(a.y, b.y, c.y, d.y) && left_turn(a, b, c) * left_turn(a, b, d) <= 0 && left_turn(c, d, a) * left_turn(c, d, b) <= 0;
};

#endif


